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1.1 Introduction 

Classical niolecular dynamics (MD) is a well established and powerful tool in 
various fields of science, e.g. chemistry, plasma physics, cluster physics and 
condensed matter physics. Objects of investigation are few-body systems and 
many-body systems as well. The broadness and level of sophistication of this 
technique is documented in many monographs and reviews, see for example [1, 
2, 3]. Here we discuss the extension of MD to quantum systems (QMD). There 
have been many attempts in this direction which differ from one another, 
depending on the type of system under consideration. One direction of QMD 
has been developed for condensed matter systems and will not discussed here, 
e.g. [4]. In this chapter we are dealing with unbound electrons as they occur 
in gases, fluids or plasmas. Here, one strategy is to replace classical point 
particles by wave packets, e.g. [4, 5, 6] which is quite successful. At the same 
time, this method struggles with problems related to the dispersion of such 
a packet and difficulties to properly describe strong electron-ion interaction 
and bound state formation. We, therefore, avoid such restrictions and consider 
a completely general alternative approach. We start discussion of quantum 
dynamics from a general consideration of quantum distribution functions. 



1.2 Quantum distribution functions 

There exists a variety of different representations of quantum mechanics in- 
cluding the so-called Wigner representation which involves a class of functions 
depending on coordinates and momenta. In the classical limit, the Wigner dis- 
tribution turns into the phase space distribution / known from classical 
statistical mechanics. In contrast to /, the Wigner function may be non- 
positive which is a consequence of the coordinate-momentum (Heisenberg) 
uncertainty. This will lead to a modification of particle trajectories which we 
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discuss below in Sec. 1.4. An important property of the distribution functions 

is that they can be used to compute the expectation value of an arbitrary 
physical observable, {A), defined by the operator A{p, q) [7], 

{A){t)= j dpdqA'^{p,q)F'^{p,q,t), 1 = J dpdqF^^ {q,p,t), (1.1) 

where A^ (p, q) is a scalar function and, for simplicity, we consider the ID 
case (generalization to higher dimensions and N particles is straightforward 
by re-defining the coordinate and momentum as vectors, q = {qi, . . . ,q7v}, 
p = {pi, . . . ,pjv}). F^ is defined via the nonequilibrium N-particle density 
operator p in coordinate representation (i.e. the density matrix), 

and A^ {p, q) is analogously defined from the coordinate represention of A. 

We now consider the time evolution of the WF under the influence of a 
general hamiltonian of the form 

N ^2 ^ 

^ = E + E ^(*) + E ^(*' (1-3) 
j=i i=i i<j 

where V{qi) and V{qi,qj) denote an external and an interaction potential, 
respectively. The equation of motion of F^ has the form [8, 7], some expla- 
nations are given in Sec. 1.4, 



dF 



w 



dt 

where the function 



+ P . VgF^ = ds F"^ {p - s, q, t) w {s, q, t) , (1.4) 



^ {s, 1't) = :^J dq'V {q - q', t) sin (^^^ (1.5) 

takes into account the non-local contribution of the potential energy in the 
quantum case. Equivalently, expanding the integral around q' = 0, Eq. (1.4) 
can be rewritten with an infinite sum of local potential terms 

OF^ , p dF^ _ ^ {h/2if'' Q2n+lpW ^ 



dt m dq ^ (2n + 1)! V ^9^"+^ ' Sp^n+i 

where (...,...) denotes the scalar product of two vectors which for an A''- 
particle system contain 3A^ components. 

If the potential does not contain terms higher than the second power of 
q, i.e. ^^|„>3 = 0, then Eq. (1.6) is simphfied and reduces to the classical 
Liouville equation for the distribution function /, i.e. 
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dt m dq dq dp' 

The Wigner funetion must satisfy a number of conditions [9], therefore, the 
initial function F^{q,p,0) cannot be chosen arbitrarily. Even if F^{q,p,t) 
satisfies the classical equation (1.7) it nevertheless describes the evolution of 
a quantum distribution because a properly chosen initial function {q,p, 0) 
contains, in general, all powers of h. In particular, the uncertainty princi- 
ple holds for average values of operators calculated with F^{q,p,0) and 
F^iq,p,t). 

One can rewrite Eq. (1.6) in the form analogous to the classical equa- 
tion (1.7) by replacing F by a new effective potential V^s defined as 

dVes dF^ _ dV dF^ d^V d^F^ ^ 
dq dp dq dp 24 dq^ dp^ 

The classical Liouville equation (1.7) can be efficiently solved with the method 

of characteristics, see e.g. Ref. [10], and this is also the basis of our QMD 
approach where an ensemble of "classical" (Wigner) trajectories is used to 
solve (numerically) the quantum Wigner-Liouville equation (1.4) which will 
be discussed in Sec. 1.4. The time-dependence of the trajectories is given by 
"classical" equations of motion 

dq^p_ dp ^ dVes{p,q,t) 
dt dt dq ' 

Of course, a direct solution of Eq. (1.9) with the definition Eq. (1.8) is only 
useful if the series is rapidly converging and there is only a small number of 
non-zero terms. 

However, there is also a principle difficulty with this approach, if the se- 
ries of terms with the potential derivatives is not converging which is the case, 
e.g., for a Coulomb potential (at zero distance). There are at least three so- 
lutions. The first is to solve the Wigner-Liouville equation by using Monte 
Carlo techniques [11]~ [14], which is discussed below in Section 1.4. The sec- 
ond is to replace the original potential on the r.h.s. of Eq. (1.8) by some model 
potential which has a finite number of nonzero derivatives as it is done e.g. 
in Ref. [15]. The third approach is to perform a suitable average of Ves, e.g. 
over a thermal ensemble of particles. This has been done both for external 
potentials and also for two particle interaction. This latter case of an effective 
quantum pair potential and its use in classical MD is disscused in the next 
section. 



1.3 Semiclassical MD simulation 



Quantum pair potentials. In order to obtain an effective pair potential which 
is regularized at zero by quantum effects (finite at zero interparticle dis- 
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tance), we consider Eq. (1.4) for 2 particles. Assuming further thermody- 
namic cquilibrimn (with a given temperature kgT = I//?), spatial homogene- 
ity and neglecting 3-particle correlations, one can solve for the two-particle 
Wigner function F^^ = F^^'^ {ruPi,r2,P2, « Ff^'^{n - r2,Pi,P2,P). This 
is now re- written as in the canonical case [7], F^^{r\ — r2,Pi,P2,0) = 
F]^'^(pi,/3)F^*^(p2,/3)e~'^^i2, which defines the desired quantum pair poten- 
tial V^. 

The first solution for was found by Kelbg in the limit of weak coupling 
[16] and has the form of Eq. (1.10) with jij 1, for details and references, cf. 
[10, 17]. The Kelbg potential (or slightly modified versions) is widely used in 
numerical simulations of dense plasmas, e.g. [18, 19, 5, 20]. It is finite at zero 
distance correctly capturing basic quantum diffraction effects which prevent 
any divergence. However, the absolute value at r = is incorrect which has 
lead to the derivation of further improved potentials, e.g. [21, 10, 17] and 
references therein. Here we use the improved Kelbg potential (IKP), 



^(n,,/3) = ^ 



+ V^Y^ ( 1 - erf 



(1.10) 



where Xij = ]rjj]/Ay, A?^ = and |J^^^ = + rrij^, which contains 
an additional free parameter jij which can obtained from a fit to the exact 
solution of the two-particle problem [17]. 

MD Simulations. We have performed extensive MD simulations of dense par- 
tially ionized hydrogen in thermodynamic equilibrium using different IKP for 
electrons with different spin projections. To properly account for the long- 
range character of the potentials, we used periodic boundary conditions with 
the standard Ewald procedure [3]. The number of electrons and protons was 
N = 200. Our MD simulations use standard Runge-Kutte or Verlet algorithms 
[3] to solve Newton's equations Eq. (1.9), where VeS is replaced by the IKP. 
Because of the temperature dependence of the IKP we applied a temperature 
scaling at every time step for all components separately (for protons and two 
sorts of electrons) to guarantee a constant temperature of all components in 
our equilibrium simulations. In each simulation the system was equilibrated 
for at least 10^ MD steps, only after this the observables have been computed. 

In Fig. 1.1 we plot the internal energy per atom as a function of temper- 
ature for two densities and compare it to path integral Monte Carlo (PIMC) 
results [17, 22]. The density is given by the Bruckner parameter = f/as, 
where f is average interparticle distance and as denotes the Bohr radius. 
For high temperatures and weak coupling, F = /{fkBT) < 1 (fully ionized 
plasma), the two simulations coincide within the limits of statistical errors. If 
we use the original Kelbg potential, at temperatures below 300,000 K (approx- 
imately 2 times the binding energy), the MD results start to strongly deviate 
from the PIMC results. In contrast the IKP fully agrees the PIMC data even 
at temperatures far below the hydrogen binding energy Ry, where the plasma 
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Fig. 1.1. Internal en- 
ergy per hydrogen atom 
at Ts = 4 and Ts = 6 
versus temperature, MD 
results are compared to 
restricted PIMC simula- 
tions [17, 22]. 



is dominated by atoms, which is a remarkable extension of 'semiclassical' MD 
into the theoretically very difficult regime of moderate coupling, moderate 
degeneracy and partial ionization. 

Interestingly, even bound states can be analyzed in our simulations by 
folowing the electron trajectories. At T < IRy, we observe an increasing 
number of electrons undergoing strong deflection (large angle scattering) on 
protons and eventually performing quasi-bound trajectories. Most of these 
electrons remain "bound" only for a few classical orbits and then leave the 
proton again. Averaged over a long time, our simulations are able to reveal 
the degree of ionization of the plasma. For temperatures below approximately 
50, OOOisr (which is close to the binding energy of hydrogen molecules), the sim- 
ulations cannot be applied. Although we clearly observe molecule formation 
(see below), there also appear clusters of several molecules which is unphysi- 
cal in the present conditions and is caused by the approximate (two-particle) 
treatment of quantum effects in the IKP. This turns out to be the reason for 
the too low energy in Fig. 1.1 at low temperatures. 

Let us now turn to a more detailed analysis of the spatial configuration of 
the particles. In Fig. 1.2 the pair distribution functions of all particle species 
with the same charge arc plotted at two densities. Consider first the case 
of T = 125,000 K (upper panels). For both densities, all functions agree 
qualitatively, showing a depletion at zero distance due to Coulomb repulsion. 
Besides, there are differences which arise from the spin properties. Electrons 
with the same spin show a "Coulomb hole" around r = which is broader 
than the one of the protons due to the Pauli principle (additional repulsion 
of electrons with the same spin projection). This trend is reversed at low 
temperatures, see the middle panel, which is due to the formation of hydrogen 
atoms and molecules. In this case, electrons (i.e., their classical trajectories) 
are "spread out" around the protons giving rise to an increased probability 
of close encounters of two electrons belonging to different atoms compared to 
two protons. 
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Fig. 1.2. Electron- 
electron and proton- 
proton pair distribu- 
tion functions for a 

correlated hydrogen 
plasma with Ts = 4 
(left row) and rs = 6 
(right row) for T = 
125, OOOA:, 61, 250if , 
and 31,250/s: (from 
top to bottom). 
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Fig. 1.3. Electron- 
proton pair distribution 
functions multiplied by 

as function of e-p 
distance at rs= 4 (top) 
and rs = 6 (bottom) and 
four temperatures. 



Now, compare electrons with parallel vs. electrons with anti-parallel spins. 
In all cases, we observe a significantly increased probability to find two elec- 
trons with opposite spin at small distances below one Bohr radius which is due 
to the missing Pauli repulsion in this case. This trend increases with lowering 
of the temperature due to increasing quantum effects. Before analyzing the 
lowest temperature in Fig. 1.2 consider the electron-proton distributions. Mul- 
tiplying these functions by gives essentially the radial probability density 
Wep(T) = r^gep(r), which is plotted in Fig. 1.3. At low temperatures this func- 
tion converges to the ground state probability density of the hydrogen atom 
^e-pir) = {r)\'4>\\s{T) influenced by the surrounding plasma. Here, lowering 
of the temperature leads towards the formation of a shoulder around lAas for 
Ts = 4 and 1.2aB for rs = 6 which is due to the formation of hydrogen atoms 
(which is confirmed by the corresponding quasi-bound electron trajectories). 
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At this temperature, the observed most probable electron distance is slightly 
larger than one ag as in the atom hydrogen ground state. Of course, classical 
MD cannot yield quantization of the bound electron motion, but it correctly 
reproduces (via averaging over the trajectories) the statistical properties of 

the atoms. 

At 62, 500K and rg = 6 (right center part of Fig. 1.2) the simulations show 
a first weak signature of molecule formation - see the maximum of the p-p 
distribution fimction aroimd r ~ 2aB and the maximum of the distribution 
function of electrons with anti-parallel spins around r = 1.5aB- Upon further 
lowering of the temperature by a factor of two (lower panel of Fig. 1.2) the 
p-p functions exhibit a clear peak very close to r = lAas " the theoretical p-p 
separation in H2. At the same time, also the e-e functions have a clear peak 
around r = O.Sas (the two electrons are concentrated between the protons). 
In contrast, in the case of parallel spins, no molecules are formed, the most 
probable electron distance is around r = 1.2aB- 

MD results for dynamic quantities. We now extend the analysis to the dy- 
namic properties of an hydrogen plasma in equilibrium which is based on the 
fluctuation-dissipation theorem. The time-dependent microscopic density of 
plasma species a is defined as p°'{r, t) = ^2^=1 ^[^ ~ ""f (*)]> '^it^i Fourier 
components /9"(k, t) = ^fj!'iCxp[ik ■ rf (f)], where ri{t) denotes the trajec- 
tory of particle "i" obtained in the simulation. We now define the three partial 
density-density time correlation functions (DDCF) between sorts a and r/ as 

A'^^ik,t) = (N^ + N,)-' (p"(k,f)p''(-k,0)) , (1.11) 

where, due to isotropy, k — fc. Here (...), denotes averaging along the trajec- 
tories by shifting the time interval and keeping the difference equal to t. Note 
also, that A°''^{k, t) ~ A^°'{k,t) for all pairs a and rj. In addition to the spin- 
rcsolvcid electron functions we can also consider the spin averaged correlation 
function A"{k,t) = A''^ {k,t) + A^^k,t). 

We have performed a series of simulation runs of equilibrium fluctuations 
in hydrogen plasmas with coupling parameters F and electron degeneracy pa- 
rameters Xe = pA'e with the elctron DeBroglie wavelength A,, = h/ \/2/iuriekBT 
ranging from zero (classical system) to one (quantum or "degenerate" system). 
The electron DDCF ior F = \ and Xe = 1 are plotted in Fig. 1.4 for four values 
of the dimensionless wavenumber, q = kf. The correlation functions (tt and 
J,t) have two characteristic behaviors - a highly damped, high-frequency part 
and a weakly damped low-frequency tail. The latter is related to the slow ionic 
motion and the first one to oscillations with frequencies close to the electron 
plasma frequency cjpi . On the other hand, the time scale of the ion motion is 
determined by the ion plasma frequency cj*; = ■\/ AirpiZfe^/rrii, their ratio be- 
ing about 43 (square root of the ion-to-electron mass ratio). The slow proton 
oscillations are clearly seen in the proton DDCF, shown in Fig. 1.5. To resolve 
the proton oscillations the whole simulation (including the electron dynamics) 
has to extend over several proton plasma periods Tp = 27r/a;L thereby resolv- 
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Fig. 1.4. Electron DDCF (1.11) mul- 
tiplied by {Nj + iVi) for r = 1 and 
Xe = 1 for four wave vectors. Upper (cen- 
ter) figure: correlation functions for par- 
allel (antiparallol) spins, bottom figure: 
spin-averaged function. 



Fig. 1.5. Proton DDCF (1.11) for 
r = 1 and Xe = 1 for five wave vec- 
tors (in units of 1/f). 



ing the fast electronic motions as well, which sets the numerical limitation of 
the calculation. 




Fig. 1.6. fon-acoustic wave 
dispersion in a dense hy- 
drogen plasma. Lines corre- 
spond to weighted linear fits 
to the MD data (symbols). 
The scatter of the data is 
due to the limited particle 
number A'^ and simulation 
time and can be systemati- 
cally reduced. Also, smaller 
q-values require larger N. 
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The temporal Fourier transform of the DDCF yields another very impor- 
tant quantity - the dynamic structure factor, Sa,ri{'jJj q), which allows one 
to analyze, e.g., the dispersion of the coupled electron and proton oscilla- 
tions. Fig. 1.6 shows dispersion results for the collective proton oscillations 
(for the electron modes, see Refs. [20, 22]) which follow from the peak posi- 
tions of Sii{LU, q). Fig. 1.6 shows the peak frequency versus wave number, i.e. 
the dispersion of longitudinal ion-acoustic waves, uj{q) = wmd • <?: where wmd 
denotes our MD result for the phase velocity. This can be compared to the 
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familiar analytical expression for an ideal two-temperature (Tg » Tj) plasma 
Vs = \/ ZikBTe/m,i. where Vg is the ion sound velocity. We observe deviations 
of about 10%, for weak degeneracy, Xe < 0.5, and aboutlO%, for large degen- 
eracy, Xe > 1) which are due to nonideality (correlations) and quantum effects, 
directly included in our simulations. For further details on this method, see 
Refs. [22, 23, 6]. 

Thus the smiclassical MD is a powerful approach to correlated quantum 

plasmas. Thermodynamic and dynamic properties are accurately computed if 
accurate quantum pair potentials, such as the IKP, are being used. 



1.4 Quantum dynamics 

Now we discuss the method of Wigner trajectories in more detail. As we have 
seen, the Wigner function W [to avoid confusion, in this section we rename 
W] in Eq. (1.2) is the Fourier transform of the non-diagonal elements 
of the density matrix which, for a pure state, is p (<? + -I, g — -1) = t/j{q + 
^,t)tp*{q— f, t), where the A^— particle wave functions satisfy the Schrodinger 
equation with an initial condition 

ih^^Hij; V (to) = (9) , (1.12) 

which contains the hamiltonian (1.3) [recall that g is a vector of dimensionality 
{Nxd}]. By taking the time derivative ofW in Eq. (1.2), substituting it in the 
l.h.s of the Schrodinger equation instead of dtp/dt and integrating by parts, 
we recover Eq. (1.4). For convenience of the further analysis, we separate the 
contribution of the classical force F{q) = —'VqV{q), compensating it in the 
function co which now replaces u), 

+ JL.^^W + F{q)-VpW^ / dsW{p-s,q,t) co{s,q,t), (1.13) 

Ot 771 J — oo 

9' *) = / ^^9'^ (1 - *) sin (^) "^sS is) . (1.14) 

In the classical limit {h 0), the r.h.s of Eq. (1.13) vanishes and we obtain 
the classical Liouville equation 

dW n 

+ IL.x^ W + Fiq)-VpW = 0. (1.15) 

ot m 

The solution of Eq. (1-15) is known and can be expressed by the Green 
function [9] 

G{p, q, t;po, qo, h) = S\p- p{t; to,po, qo)]S[q - q{t; to,po, go)], 

where p{t) and q{T) are the phase space trajectories (of all particles), which 
are the solutions of Hamilton's equations together with the initial conditions 

at T = to = 0, 
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dq/dr = p{T)/m; g(0) = qo, 

dp/dT = F{q{T)y, p{0)=po. (1.16) 

Using the Green function, the time-dependent solution of the classical Liou- 
ville equation takes the form 

W {p,q,t) = j dpodqoG{p,q,t;po,qo,0)Wo{po,qo)- (1-17) 

With this result, it is now possible to construct a solution also for the quan- 
tum case. To this end we note that it is straightforward to convert Eq. (1.13) 
into an integral equation 

W {p, q,t)= j dpodqo G{p, q, t;po, qo, 0) Wo(Po, ^o) + (1-18) 

/ dti I dpidqiG{p,q,t;pi,qi,ti) / dsiLo{si,qi,ti)W{pi - Si,qi,ti), 

Jo J J —oo 

which is exact and can be solved efficiently by iterations [11, 10]. First, wc 
can express the function W{p\_ — si,q\,t\) written on the r.h.s. of Eq. (1.18), 
formally, using the same Eq. (1.18) but written for another set of variables, 
i.e. {p, g, t} {pi - si,qi,ti] and {pi,qi, si,ti} {p2, ^2, S2, i2}. Second, 
substitution of the obtained result again in Eq. (1.18) gives 

W {p, q, t) = (p, q, t) + {p, q, t) + dti J dl G{p, q,t;l,h) x 

rti p roc 

/ dt2 d2G{pi-Si,qi,ti;2,t2) ds2U){s2,q2,t2)W{p2 - S2,q2,t2), 

Jo J J -oo 

(1.19) 

where we have introduced the short notations n = qn,Pn, dn = dqndpn and 
W^°\p,q,t)= J dOG{p,q,t;0,0)Wo{0), (1.20) 
W^^^p,q,t)= dh dlG{p,q,t;l,h) dsiuj{si,quti) 

•/ 'J —oo J — oo 

X J dOG{p,-suqi,ti;0,0)Wo{0). (1.21) 

W^^^ {p, q, t) (as it follows from the Green function G{p, q, t; po,qo, 0)) describes 
the propagation of the Wigner function along the classical characteristics, i.e., 
the solutions of Hamilton's equations (1.16) in the time interval [0,t]. It is 
worth mentioning, that this first term describes both classical and quantum 
effects, due to the fact that the initial Wigner function Wo(poj So), in general, 
contains all powers of Planck's constant h contained in the initial state wave 
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functions. The second and third terms on the r.h.s. of Eq. (1.19) describe 

additional quantum corrections to the time evolution of W{p, q, t). 

Let us consider the term W'^^^ (p, q, t) in more detail. It was first proposed in 
Ref. [11] and demonstrated in Refs. [12]~[14] that the multiple integral (1.21) 
can be calculated stochastically by Monte Carlo techniques. For this wc need 
to generate an ensemble of trajectories in phase space. To each trajectory we 
ascribe a specific weight, which gives its contribution to (1.21). For example, 
let us consider a trajectory which starts at point {pQ,qQ^T = 0}. This trajec- 
tory acquires a weight equal to the value Wo(Po,9o)- Up to the time t = tx 
the trajectory is defined by the Green function G{p\ — s\,q\,t\;po, qo, 0). At 
T = ti, as it follows from Eq. (1.21), the weight of this trajectory must be 
multiplied by the factor a;(si, gi, fi), and simultaneously a perturbation in 
momentum takes place: {pi — si) Pi- As a result the trajectory becomes 
discontinuous in momentum space (but continuous in the coordinate space). 
This is, obviously, a manifestation of the Heisenberg uncertainty of coordi- 
nates and momenta. Now, the trajectory consists of two parts - two classical 
trajectories which arc the solutions of Eq. (1.16), which arc separated, at 
T = ti, by a "momentum jump" of magnitude si. What about the value si of 
the jump and the time instance t\l Both appear under integrals with a cer- 
tain probability. To sample this probability adequately, a statistical ensemble 
of trajectories should be generated, further the instant of time tx must be 
chosen randomly in the interval [0, t], and the momentum jump s\ randomly 
in the interval (— oo,-|-oo). Finally, also different starting points {poj^o} of 
trajectories at r = must be considered (due to the integration / dpodqo). 
Considering a sufficiently large number of trajectories of such type we can ac- 
curately calculate W^^'' (p, q, t) - the first correction to the classical evolution 
of the quantum distribution function W^'^^p, q,t). 

Let us now take into account the third term in Eq. (1.19). We now sub- 
stitute, instead of W{p2 — 52,92,^2), its integral representation, using (1.18). 
As a result we get, for this term, 

W^'^\p,q,t)^ / dh / dlG{p,q,t;lM) / dsi a;(si, 91, ti) x 

rti p PCX} 

/ dt2 / d2G{pi - si,qi,ti;2,t2) / ds2 w(s2, ^2, ^2) x 

Jo J J-00 

dOGip2 - S2, q2M; 0, 0) Wo(0). (1.22) 



If we apply the stochastic interpretation of the integrals, as we did above for 
W^^'> {p,q,t), this term can be analogously calculated using an ensemble of 
classical trajectories with two momentum jumps taking place at time instants 
T = t\ and T = t2, and with a weight function multiplied by the factors 
w{s\,q\,t\) and ui{s2,q2,t2), respectively. 

Applying the above procedure several times, we can get the higher order 
correction terms. As a result, W{p, q, t) will be expressed as an iteration series, 
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with each term of the series representing a contribution of trajectories of 

a definite topological type - with one, two, three, etc. momentum jumps. 
In Fig. 1.7 we show an example of trajectories contributing to the terms 
and 




Fig. 1.7. Illustration of the iteration series. Three types of trajectories are shown: 
without (top curve), with one (middle) and two (lower) momentum jumps. 

As was noted in Section 1.2 the Wigncr fimction allows us to compute 
quantum-mechanical expectation values of an aribtrary one-particle operator 
A. Using the idea of the iteration series (1.19), we obtain an iteration series 
also for the expectation value, 

{A)it) = J dpdqA{p,q)W{p,q,t) = {A)(°\t) + {A)^'\t) + ..., (1.23) 

where different terms correspond to different terms in the series for W. The 
series (123) maybe computed much more efficiently than the one for W since 
the result does not depend on coordinates and momenta anymore. 

Certainly, in the iteration series it is possible to take into account only a 
finite number of terms and contributions of a limited number of trajectories. 
Interestingly, it is not necessary to compute the individual terms iteratively. 
Instead, all (relevant) terms can be calculated simultaneously using the ba- 
sic concepts of Monte Carlo methods, e.g. [24]. An important task of the 
Monte Carlo procedure will be to generate stochastically mostly the trajecto- 
ries which give the dominant contribution to the result, for details, see [10]. 

1.5 Time correlation functions in the canonical ensemble 

So far we have considered the dynamics of pure states where the density matrix 

p is defined by a single wavefunction ip. However, at finite temperature p is, 
in general, defined by an incoherent superposition of wave functions (mixed 
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states). Here we consider a canonical ensemble as the most common one. Time 
correlation functions CpAit) = {F{0)A{t)) arc among the most important 
quantities in Statistical physics which describe transport properties, such as 
diffusion, dielectric properties, chemical reaction rates, equilibrium and non- 
equilibrium optical properties, etc. An example has already been considered 
in Sec. 1.3 - the density-density auto-correlation function (1.11). Here, we 
use a more general expression for the quantum correlation function of two 
quantities A and F given by the operators F and A. In the canonical ensemble 
the averaging is performed by a trace with the canonical density operator 
^EQ _ z~^e~^^ , with (3 — l/ksT, and the correlation function has the 
form [25] 

CFA{t) = ^Ti{Fe'"''f< ie"*-^*"), (1.24) 

where Fi is the Hamiltonian (1.3), is a complex time argument (it "absorbs" 
P^^)i = t — i/3/2, Z = Tr/9^'3 jg the partition function, and we use h= 1. 

The time correlation function can now be computed by first, writing 
Eq. (1.24) in coordinate representation and then transforming to the Wigner 
picture using the Weyl representation of F and A, 

CpAit) = 

Z-^ J dqidq2dqsdq^ {qi\F\q2) {q2W"*^qi) {qi\A\qi) {qi\e-'"'^\qi'^ 
dpidqidp2dq2 F{pi,qi) A{p2, 92) W{pi, qi;p2, q2] t; (3), (1.25) 



where W{pi,qi;p2,q2]t: 0) is now a generalization of the Wigner function 
which is defined as double Fourier transformation of the product of two non- 
diagonal matrix elements of the density operator 

Wipuqi;P2,q2;t;P) = -^^^^ J d^id^2e'P'^' e*f^«^ x 



6 



^C2\/ 6 



+ (1.26) 



Calculating the partial time derivatives of the function W using the usual 
technique for matrix elements, it can be shown that the function W satisfies 
a system of two Wigner-Liouville equations [12, 13] 



at m 
dW m 

_ + ^ . Vg^W + F{q2) ■ V^^W = I2, (1.27) 



where on the r.h.s. we have two "coUison" integrals 
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/i = / dsiW{pi- Si,qi;p2,q2;t;(3) oj{si,qi,t) 

J —OO 

/oo 
ds2W {pi,qi;p2 - S2,q2;t; l3) uj{s2,q2,t), (1-28) 
-oo 

and the function oj (s, q, t) is defined in the same way as in the microcanonical 
ensemble, see Eq. (1.14). 

1.5.1 Initial conditions for the Wigner-Liouville equation 

Using Eq. (1.26) at i = 0, wc find that the initial value of the Wigner function 
is given by the integral (again I = qi,pi') 



Wo(l;2;0;/3) 



^1 



Z(27r) 

-0H/2 



' DH/2 



C12 + ^)U2-^ 



gi + |>. (1.29) 



Let us now use the groiip property of the density operator p and the high 
temperature approximation for the matrix elements of {q'\p\q) [26], 



-I3H _ 



M 



e 2M 



H 



rK 



(1.30) 



As a result, we obtain 



Wo(l;2;0;/3) 



1 



/ 
/ 



Z(27rfi) 



2Nd 



dq[ . . . dq'M dq'l . . . dq'li e ".=2 m=i x 



rK 



91 + 1- ){qi 

92 + y ) {q2 



2 
2 



rK 



-J^K 



(1.31) 



where 



Here we have assumed that M ^ 1, and — denotes the thermal de 

Broglie wave length corresponding to the inverse temperature ^fy. A direct 
calculation of the last two factors in Eq. (1.31) gives 



/ 



rK 



6 



91 



ei 



g 2M 



K 



rK 



91/ , 



(1.32) 
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where 

4>{P\<lM,1i) = {^Xm)""'^ exp I — 1 . (1.33) 

The final result for the Wigner function at t = can be written as 

W^(l;2;0;/3)« 

j dq[...dq'j,,dq'l ...dq';,^{l;2;q[ 
(l>{P2;q'M,<li)HPi'':QM,Qi), 

where 

M + l M 

^{pi,qi]P2,q2]qi---qM-^qi ■■■qll^P) = Z -=1 -=1 . 

Here we have introduced the notation {^q = gii^o = ^2} and {g^+i = 
Q2;<?M+i — The Fig. 1.8 is an illustration of the simulation idea. Two 
closed loops with the set of points show the path integral representation of 
the density matricies in Eq. (1.31). The left chain of points, i.e. 
{qi,qi, . . . ,q'fj,q2,qiT . ■ , q'l^} characterizes the "path" of a single quantum 
particle. The chain has two special points {pi,qi) and {p2,q2)- As it follows 
from Eqs. (1.26), (1.27), these points are the original points for the Wigner 
function (the additional arguments arise from the path integral representa- 
tion). As we show in the next section, we can consider these points as starting 
points for two dynamical trajectories propagating forward and backward in 
time, i.e. f — > t+ and t ^ t~ . The Hamilton equations for the trajectories are 
defined in the next section. 



••■9m;9i •••9m;0;/3) X 

(1.34) 



1.5.2 Integral equations 

The solution scheme follows the one explained before. The only difference is 
that we now have to propagate two trajectories instead of one, 

^ = -^FMt)]; P2iO)=pi (1.35) 

the first (second) propogating forward (backward) Let us substitute expres- 
sions for F[qi{T)],pi{T), F[q2{T)] and ^2(7') from (1.35) into Eqs. (1.27) and 
subtract the second equation from the first. As a result, in the l.h.s. we will 
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q, 




/' t" 



\ t" 



Fig. 1.8. Two closed loops illustrating the path integral representation of two elec- 
trons in the density matricies in Eq. (1.31), see also Ref. [26]. Two special points 
(j>i,qi) and (^2,92)- are starting points for two dynamical trajectories propagating 
forward and backward in time. 

obtain a full differential of the Wigner function. After multiplication by the 
factor {1/2} and integration over time, the integral equation for the Wigner 
function takes the form 



W {pi,qi;p2,q2;t;(3) = 

I dpldqldp'ldql G(pi,(7i,p2,g2,t;p?,g?,P^,92°,0) M/(p;,g?;p°,gO;0;/3) + 



w\icvc "d{s,q\:ri,q\;T) = ^{ijj[s,q\)5{'q) — uj{ri,q\)5{s)}. The dynamical Green 
function G is defined as 



G(pi,gi,P2,g2,i;p?,g?,P2,92,0) = 5{pi- pi{T]p\,q^,())] x 

5{qi - qi{T;pl, g?,0)} • 5{p2 - P2{t;P%, ql, 0)} • 5{q2 - q^T^p^ ql, 0)} (1.37) 



Let us denote the first term on the r.h.s. of Eq. (1.36) as W^^^ {p\,q\]p2, (72; t; (5). 
This term represents the Wigner function of the initial state propagating along 

classical trajectories (characteristics), solutions of Eqs. (1.35). Using the same 
approach which was applied in the microcanonical ensemble, we obtain ex- 
pressions for 



and represent (pi, (j'i;p2, (72; i; /3) as an iteration series. In this case, we 
can calculate this also with an ensemble of trajectories using the Quantum- 
Dynamics- Monte-Carlo approach which was described in [27]. As a result the 




(1.36) 



W^^\pr, qi;P2, q2; t; /3), W^^^ipuqi;P2, 92; t; l3),... 
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expression for the time correlation function (1.25) can be rewritten as 



where (. . . | . . .) denotes the integral in the (now 2A''-particle) phase space 
{Pi,qi,P2,q2}, and (/)(P) = F{pi,qi) A{p2,q2). 

An illustrative example for the calculations of the time correlation func- 
tions CpA is the momentum-momentum autocorrelation functions Cpp{t) for 
a one-dimensional system of interacting electrons in an array of fixed random 
scatterers at finite temperature [27]. This system is of high interest because 
at zero temperature it shows Anderson localization if eIectron-ele(;tron (e-e) 
interaction is neglected. It has been a long standing of the question what the 
effect of e-e interaction on localization will be. The present method [27] is 
well suited to answer this question, and it could be shown that Coulomb e-e 
interaction enhances the mobility of localized electrons [27, 10]. 

1.6 Discussion 

We have presented a general idea how to extend the powerful method of molec- 
ular dynamics to quantum systems. First, we discussed "semiclassical MD", 
i.e., classical MD with accurate quantum pair potentials. This method is very 
efficient and allows to compute thermodynamic properties of partially ion- 
ized plasmas for temperatures above the molecule binding energy (i.e. as long 
as three and four particle correlations can be neglected). Further, frequency 
dependent quantities (e.g., plasmon spectrum) are computed correctly for 
w < Wpi. Further progress is possible if more general quantum potentials are 
derived. 

In the second part, we considered ways to rigorously solve the quantum 
Wigner-Liouville equation for the N-particle Wigner function. Results were 
derived for both, a pure quantum state and a mixed state (canonical ensem- 
ble). Although this method is by now well formulated, it is still very CPU 
time costly, so that practical applications are only starting to emerge. Yet, 
we expect that, due to its first principle character, Wigner function QMD 
will become increasingly important for a large variety of complex many-body 
problems. 

This work is supported by the Deutsche Forschungsgmeinschaft via SFB- 
TR 24 and in part by Award No. Y2-P-11-02 of the U.S. Civihan Research 
and Development Foundation for the Independent States of the Former Soviet 
Union (CRDF) and of Ministry of Education and Science of Russian Feder- 
ation, and RF President Grant NS-3683.2006.2 for governmental support of 
leading scientific schools. 




dpidq^ dp2dq2 F{pi, qi) A{p2, 92) W{pi,qi;p2, q2;t; /?) = 



00 



(0(P)|^^(°)(P;/3)) +£(0(P)|V^«(P;/3)) 



(1.38) 
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